Deterministic endless collective evolvement in active nematics 
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We propose a simple deterministic dynamic equation and reveal the mechanism of large-scale 
endless evolvement of spatial density inhomogeneity in active nematic. We determine the phase 
regions analytically. The interplay of density, magnitude of nematic order, and nematic director 
is crucial for the long-wave-length instability and the emergence of seemingly fluctuated collective 
motions. Ordered nematic domains can absorb particles, grow and divide endlessly. The present 
finding extends our understanding of the large-scale and seemingly fluctuated organization in active 
fluids. 



Assemblies of active particles, which may be a 
physical abstraction of running animals [lL hying 
birds 0, swimming bacteria migrating cells[4[ or even 
cytoskeletonjlj, have been served as a new building block 
for physicists over the last decade or so to understand 
the common collective behavior in these non-equilibrium 
systemsQ. Active nematic, a recently proposed concept 
for a kind of apolar active particles, is formed by driven 
rod-like particles with head-tail symmetry through the 
randomly driving force along the rod orientation axis at 
the single-particle level[7(a)-7(f)]. The symmetry of the 
system is not broken by applying such micro-driven forces 
until spontaneous symmetry breaking occurs. Recently, 
simulations and experiments in active nematic system 
show well-organized collective motions with system-sized 
fluctuation[7(b)-7(e)]. For example, splitting and merg- 
ing of large-scale self-organized structures are exemplified 
in the simulation on active nematics [7(c)]. Experiments 
also show large-scale collective swarming and swirling in 
driven granular rods monolayer [7(d), 7(e)]. These obser- 
vations lead us to think about the nature of such seem- 
ingly fluctuated collective motions. It is currently un- 
clear whether these large-scale collective motions arise in 
a deterministic manner or as a result of noises applied 
upon the system. Moreover, are they genuinely restless 
on large scale and evolving without end? These impor- 
tant aspects are still not well addressed in previous stud- 
ies. 

In the present study, we start from a deterministic 
equation to study the mechanism of restless collective 
evolution in active nematics. We reveal a new phase that 
is characterized by the unattainability of stable steady 
state. We first identify that, if the steady state can be 
reached, the system investigated here favors spatially ho- 
mogeneous state. On the other hand, by taking account 
of the interplay between particle density and local ne- 
matic order (i.e., magnitude and orientation), the linear 
stability analysis shows that homogeneous nematic state 
can be unstable to fluctuations of small wave number. 
Therefore, the system enters into a chaotic phase region 
with no stable steady state. Large-scale spatial inhomo- 
geneity of density and nematic order is developed as a 
result of long-wavelength instability. The spatial inho- 



mogeneity in turn changes the direction of the nematic 
director, leading to a non-ending evolvement of the sys- 
tem. Numerical flux analysis shows that the particle-rich 
nematic domains are surrounded by particles fluxes, and 
evolve via absorbing particles from low-density isotropic 
medium, growing, and extending itself and breaking into 
small pieces. More importantly, all these seemingly fluc- 
tuated collective motions giving rise to giant number fluc- 
tuations are governed by a deterministic equation which 
is essentially free of noises. 

We notice that one salient feature of simulation rules 
for active nematics by Chate et al.[7(c)] is that particle 
rotations are governed through inter-particle nematic in- 
teraction while spatial translational movements are free 
of such interactions. Experimentally particles arc driven 
along their long axis, inducing strong longitudinal diffu- 
sion, and they can thrust into the media with the supply 
of kinetic energy [7(e)]. A simple diffusion equation which 
follows these observations can be written as (see ||): 



d t f(r, u, t) = V[D||UuV + D ± (I- uu)V]/(r, u) 
+ &[Dr&f{T, u) + D r .%w(r, u)/(r, u)] , 



(1) 



where Du and D± are the parallel and perpendicular 
components of the translational diffusion constants. D r 
is the rotational diffusion constant, and the rotational 
operator is defined by = u x 9 U 0. /(r, u, t) is the 
particle number distribution function where the spatial 
coordinate r and the unit vector u denote the center- 
of-mass position and long-axis direction of particles, re- 
spectively. w(r, u) is a self-consistent interacting poten- 
tial which has ±u-symmetry. In two-dimensional(2-D) 
case, the most common form of such interacting po- 
tential is the excluded-volume-like interaction w(r, u) = 
I 2 J du'|u x u'|/(r, u'), where I is the particle length. 

The diffusion equation Eq.([T|) for active nematics 
satisfies particle number conservation with the spa- 
tial translational current J s (r, u) = — _DiiuuV/(r, u) — 
D±(I — uu)V/(r,u) and the local rotational current 
J r (r,u) = -D r Mf (r,u) - D r ^w(r,u)f(r,u), which 
are independent of each other. The translational cur- 
rent is purely diffusive. In the Fourier space as de- 
fined by /(r, u) = J dr/(k, u)e~ lkr , the spatial fluc- 
tuation modes are governed by diffusive decaying term 
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FIG. 1: Three phases are separated by the blue curves. Polar 
plots (a)-(f) indicate the instability regime for k and 9. The 
locations (Do, S) of their origins are used as the parameters 
to produce these plots, (a) \D \ = 0, 5 = 0.8, (b) \D \ = 2/3, 
5 = 0.4, (c) \D \ = 0, 5 = 0.4, (d) |D | = 1/3, r5 = 0.01, 
(e) | A) | = 2/3, 6 = 0.01, and (f) |A>| = 0, 5 = 0.01. The 
black and red branches in polar plots represent the cases of 
positive and minus Do, respectively. The inset shows the 
corresponding instability modes D^ 1 gp |e=o for these polar 
plots. 



— fc 2 (Z?||Cos 2 </9 + Dj_sin 2 ( / j)/(k, u), where tp is the angle 
k makes with u. For positive Dm and D±, all the spa- 
tial fluctuation modes decay except for k — which in- 
dicates total particle number conservation. Therefore, 
as governed by these decaying diffusive modes, the spa- 
tial term suggests that only the spatially homogeneous 
state will probably be the stable steady state if there 
are no particle sources and sinks in the system and at 
the boundaries. However, when D\\ ^ D±, we will show 
that the spatially homogeneous state may become unsta- 
ble to fluctuations in the nematic state. Consequently, 
the system is deprived of all possible stable steady state 
and becomes restless and evolves endlessly, similar to the 
deterministic nonperiodic flow found in turbulence [loj. 

We first examine the spatially homogeneous dynamic 
equation derived from Eq.{T]) by neglecting the spatial 
derivatives: d t f(u,t) = 3g[D r 3!f(u) + D r Mw(u)f(u)}. 
Determined by the integration kernal of the self- 
consistent interacting potential w(u), f(u,t) has ±u- 
symmetry. In spatially homogeneous nematic state, we 
assume that the nematic director no is in the x-axis, 
and thus the distribution function can be expanded as 
f(u,t) = (2n)~ 1 pJ2^ =0 a n (t)cos2n(j), where p is the par- 
ticle number density, (f> is the angle of the unit vec- 
tor u, and n = 0, 1,2 •■•oo. The dynamic equation 
of the coefficient a n (t) is given by (4D r ) -1 <9 t a„(t) = 

2 , pl 2 n 2 a n , ,2 y^°o nma m (a, , -a\ n+ „ A ) , 

n an+ (4n i_ 1)T +f« Z^m=i (4m a -i)7r , wneie 

ao(t) = 1 and ai(t) = 2S(t), since the number den- 
sity p — J duf(u,t) and the nematic order parameter 
S(t) = Jducos20/(u, t)/p. By setting a n = for n > 3, 



the truncated dynamic equation for the nematic order 
parameter can be written as: dt 4 ^^ = (p — 1)S — J^fdrjjyj 
where p = pj p* is the rescaled number density, and the 
critical density p* = 3tt/21 2 , beyond which the system 
enters into a spatially homogenous nematic state. 

Next, we examine the linear stability of such spatially 
homogeneous nematic state above p* , by expanding the 
number distribution function /(r, u) = (27r) _1 p(r)[l + 
A{u a up — 5 a /3/2)Q a p(r)] with the inclusion of the align- 
ment tensor Q a p(r) = S(r)(n a np(r) — 5 a p/2) where n(r) 
is the unit vector of the nematic director. We assume 
that nematic director is along the x axis of the system. 
Small fluctuations of density and nematic director near 
the ordered nematic state are given by 5p(r) = p(r) — po 
and S5n y (r) = Qxy(*), respectively. Here, po = po/p* 
is the reduced bulk particle density, and Sn y (r) is the 
y-component of the deviated nematic director. Noting 
that n (r) = x and |n| = 1, Sn y (r) is the only possible 
small fluctuation of nematic director n(r). The resulting 
hydrodynamic equations can be obtained from Eq. |T]), 
yielding 

d t 5p = ^d 2 Jp+^(d 2 x -d 2 )(pS) 



p Sd t 5n y 
dt[pS(r)} 



r, -a-r ' 2 y " x "V 

+2D n p Sd x d y 5n y , 
" d x d y Sp + ^f-poSdl5n y , 
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(dl-d 2 y )p+^dl(pS) 
-AD r ~p{~p-l)S-D r ^A,., 



(2) 
(3) 

(4) 



where D p = D\\ + D± and D n = Dn — D±. It is 
easy to see that the homogeneous state is stable to fluc- 
tuations of coupled nematic director and density field. 
Here we consider the stability of the modes that cou- 
ple fluctuations of density Sp(r) = p(r) — po and mag- 
nitude of nematic order SS p (r) = p(r)S(r) — PqSq with 
5n y = around the homogeneous state (po,So), where 
So = \/4(5 — po)(po — 1)/3/5q. The mode of fluctuations 
in Fourier components with wave vector k, defined by 
5p(r) = J dkp k e- <k r and 6S p (r) = J dkSp ke - 4kr , is 
governed by 



dt 



S p k 



D p k 2 D n cos20fc 2 ' 



D n cos20k 2 /2 D p k 2 
-8D r p S +16D r r5 



S p w 



(5) 



where 8 = p — 1, 9 is the angle between wave vector k 
and nematic director no. The eigenvalues of the coeffi- 
cient matrix in Eq.© are given by D~ X X^ s = — (8r5 + 



k 2 )/2 ± yj Dq k 4 cos 2 2<9/8 - ActDqk 2 cos26» + 16<5 2 , where 
(7 = J (4 — 5)5/3, the rescaled coefficient Dq = D n /D p , 
and the wave number k = <J 'D p /D r k. The real part 
of \J S is always negative, representing stable decaying 



mode. However, the mode \~~ s becomes positive when 
32(A)<7 cos26» + S)k 2 + (2 - D\ cos 2 26»)k 4 < 0. The coef- 
ficient of k 4 is always positive since \Do\ < 1, signifying 
that for large enough wave numbers, the fluctuations are 
always stable. For small wave numbers which describe 
large-scale fluctuations, the stability is controlled by the 
coefficient of k 2 . Thus when (D aco(Q8 + S) < 0, the 
system becomes unstable on large scale. 

The phase map for (Do, 5) is given in Fig. 1. Be- 
tween the isotropic and linearly stable nematic phases, 
there is a region where spatially homogeneous nematic 
state is unstable. It is denoted as a 'phase with no sta- 
ble state' to emphasize that the only possible form of 
steady-state solution is unachievable there. For different 
(Do, 6) within the 'no stable state' region, the instabil- 
ity mode structures (k, 9) are given by the polar plots 
whose central positions represent (Do, 6). Here, each po- 
lar plot is composed of horizontal (black) and vertical 
(red) branches enclosing unstable fluctuation modes, cor- 
responding to Do < and Do > 0, revealing that spatial 
inhomogeneities are developed parallel and perpendicu- 
lar to the nematic director, respectively The maximum 
values K m of k for the instability regimes are always in 
the directions 9 = ir/2,3ir/2 for Do > and 9 = 0, 7r 
for Do < 0, respectively. Generally speaking, K m be- 
comes larger when (Do, 5) is far from the phase bound- 
ary. The inset of Fig. 1 shows the value of Z?~ 1 A^" s , 
where for small k, D^X- s > corresponds to the long- 
wavelength instability and the onset of large-scale spatial 
inhomogeneity. 

What will happen in the phase region where there is 
no stable steady state? And how the system evolves in 
time? To answer these questions we directly integrate 
Eq. ([T]) numerically in this region using alternative im- 
plicit algorithm(see [§]). Starting from an isotropic and 
spatial-homogeneous initial condition, local ordered ne- 
matic domains form at the beginning, accompanied with 
quick development of density inhomogeneity. Further 
coarsening of these structures leads to the coexistence 
of particle-enriched nematic domains and particle-poor 
isotropic region where p < p* (Fig. 2a). However, such 
a large-scale spatially inhomogencous structure is unsta- 
ble, and it will evolve and become fragmented as shown 
in Fig. 2b. The fragmented structure will again coa- 
lesce and similar process will repeat aperiodically and 
endlessly (see Q Ml.mov). In Fig. 2c-2k, we show how 
a particle-rich nematic branch breaks into pieces and re- 
unites into a structure with new morphology. 

How does the fragmentation process occur? In Fig. 
3a-c, we take a close look at the process that a nematic 
band breaks up(for a more continuous process, see [|| 
M2.mov). In Fig. 3a, after the spontaneous formation 
of a nematic band, initially, it is shown that the nematic 
director in the high-density ordered region is almost par- 
allel to the density stripe boundary. In this case, the 




FIG. 2: Density and nematic order profiles are plotted. The 
color scale shows the local relative rescaled density value 
5(r) — p(r) — 1. The length and angle of white segments 
show the magnitude and direction of nematic order, respec- 
tively. The dynamic parameters are D r = 2, D± = 0.4, 
and Du = 2.4. The reduced instability dynamic parameter 
Do = 5/7 and S = 0.01. The system size is 300 x 300 with 
the particle length 1 = 1. Periodic boundary condition is im- 
plemented. Discrete time step A t = 0.018 and spatial steps 
A x — A y = A r = 3. (a)-(b) The snapshots are taken at 
times 1.7 x 10 6 A t and 2.3 x 10 6 A t . (c)-(k) A close look of 
the breaking and coalescing processes of a nematic domain, 
at times 1.6 x 10 6 , 1.7 x 10 6 , 1.71 x 10 6 , 1.72 x 10 6 , 1.73 x 10 6 , 
1.74 x 10 6 , 1.75 x 10 6 , 1.76 x 10 6 and 1.8 x 10 6 in unit A*. 



nematic director is along the x-axis. For D n > 0, from 
our previous stability analysis, the term — ^D n d 2 (pS) 
in Eq.((3]) is directly responsible for the development of 
such density inhomogeneity. Now, we are interested if 
such aligned director field is stable to small fluctuations 
<5nj_(r) = n(r) — no = (0, Sh± y ). To linear order, the dy- 
namic equation for 6n± y (r) can be obtained from Eq.([T]) 



as pSd t Sn± y = \D p [pSd 2 8r 



.jdypS], where we 



have assumed that there is no spatial variation of pS 
along x-axis. The spatial variation of pS along the y-axis 
is significant since density inhomogeneity is developed in 
that direction. Near stripe boundaries, we always have 
d 2 pS > 0, which makes the fluctuations Sn± y unstable. 
Such instability will induce the change of the nematic 
orientation, and this explains why the nematic directors 
in Fig. 2 and Fig. 3b arc most likely to be oblique to 
the density profile boundaries. When the nematic direc- 
tors become oblique to the boundary, as shown in Fig. 
3b, there is a leakage of particles from the high-density 
region. As the particle density in the stripe falls into the 
'no stable state' region as shown in Fig. 1, the spatial 
instability takes place again. This leads to a fragmen- 
tation event, as shown in Fig. 3c. It is shown that a 
density crevice forms parallel with the nematic director 
as spatial instability requires. In Fig. 3d, we show a 
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FIG. 3: (a)-(c) Three typical steps that an ordered stripe 
breaks, at times 30000A t , 60000A t and 63000A t in sequent, 
(d) Twisted spindle shaped structure formed after it breaks off 
from a larger ordered structure. The color scale shows the lo- 
cal relative rescaled density value S(r) = p(r) — 1. The length 
and angle of white segments show the magnitude and direc- 
tion of nematic order, respectively. The black arrow shows 
the direction and strength of particle fluxes, (e)-(f) Simula- 
tion result shows density inhomogeneity and similar restless 
evolvement in active nematics for 9.4 x 10 4 and 1.1 x 10 5 
simulation sweeps. 



under a self-organized rachct potential|ll|. When the 
crevice forms in Fig. 3c, particles flow into the low den- 
sity region as guided by the nematic directors, accom- 
panied with the growth of nematicly ordered tips near 
the crevice. In Fig. 3d, as indicated by the density cur- 
rents, the twisted-spindle shaped nematic region absorbs 
particles from the medium on both sides and generate 
outward flux on both tips. In this way it can grow and 
extend itself quickly into the low-density medium. With 
the presence of these density fluxes around the ordered 
structure, it behaves like a creature absorbing, growing, 
dividing and dissipating into isotropic medium endlessly. 

In summary, the dynamic equation abstracted from 
previous simulations and experiments suggests a nemat- 
icly ordered phase with no stable steady state. Thus the 
system must evolve endlessly. We reveal the statistical 
mechanism that governs the large-scale and seemingly 
fluctuated collective evolution. We show that, as a result 
of long- wavelength instability, density and order inhomo- 
geneity develops as guided by nematic director field. The 
spatial inhomogeneity in turn changes the directions of 
local nematic directors. The changed nematic directors 
further guide the fragmentation events, leading to endless 
evolution of the system. Finally, it would be interesting 
to extend our analysis to other active fluids. 

This work was supported by the National Natural Sci- 
ence Foundation of China (No. 10974080). 



twisted-spindle shaped high-density region with local ne- 
matic order, which is commonly formed after it breaks 
off from a larger ordered structure in Fig. 2e. 

We also perform simulations to examine the stability 
of homogeneous nematic state on the basis of Eq. (JXJ) (see 
Q). Fig. 3e shows the formation of a particle-enriched 
nematic band. We find that such band is also unstable 
and undergoes similar process (see [|[ M3.mov ), where 
the nematic director changes its direction in time and the 
fragmentation event takes place afterward (Fig. 3f). 

It is worth to notice that all these highly dynamic 
structures are surrounded by particle fluxes around 
the density profile boundaries. The density currents 
arc defined as J(r) = ( J x (r), J y (r)), with J Q (r) = 
~^D p d a p(r) — D n dp[p(r)Q a p(r)] where the first term 
is just ordinary diffusive current and the second term 
is the current generated by coupling nematic directors. 
As we show in Fig. 3a, inward currents, generated by 
(0,^D n d y (pS)) which is included in the second term 
of J(r), are directly responsible for the development of 
density inhomogeneity. As the time evolves, in Fig. 
3b, along the two sides of stripe boundaries the sys- 
tem generates anti-parallel currents which also originate 
from — D n dp[p(r)Q a p(r)], and the particles seem to move 
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